#	 (C) Copyright 2007 Jason Kim Chong Polak
#	 You may reach me by email at automorphism (at) gmail (dot) com
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
#############################################################################
#Paramters
#Tweak these values according to experiment

# Analysis Parameters
autorun = TRUE
show_mouse_number = TRUE
mouse_number_format_remove = "....-..-..-"
ftype = ","

# Data Parameters
meta_data_length = 17
start_time = 0
total_time = 1500

# Rearing
r_minimum_surface_area_diff = 38
r_samples_used = 2

# Time In Center
c_distance_from_edge = 3

#############################################################################
####
#Primitive Functions
####

distance <- function(a,b)
{

#returns distance between two points. This is an internal function and is only used in other functions.
	dist = 0;
	if (is.na(a[2]) || is.na(a[1]) || is.na(b[2]) || is.na(b[1]))
	{
		dist=0;
	}
	else
	{
		dist=sqrt((a[2]-b[2])^2+(a[1]-b[1])^2);
	}
	return(dist)
}

####
#Statistics
####

maxSpeed <- function(x,y,timepoints)
{

#mainly for testing integrity of data. Most of the data has a few errors here and there and thus this function is not very reliable.
	
	max = 0;
	newS = 0;
	for (i in 1:(length(x)-1))
	{
		p1 = c(x[i],y[i]);
		p2 = c(x[i+1],y[i+1]);
		newS = distance(p1,p2)/(timepoints[i+1] - timepoints[i]);
		if (newS > max)
		{
			max = newS;
		}
	}
	
	return(max);
}

verticalExposure <- function(sa)
{

#this function simply returns the average surface area of the mouse. The more times the mouse rears or hides, the less this average will be. However, there are some real problems with the data and using this function because of the position of the camera. However, rearing will still affect this average, but the overall absolute average will definitely not be the same across all mice because of the assymetry in the boxes.

	return(sum(sa,na.rm=T)/length(sa));
}

rearing <- function(sa,mca,samples)
{
	freq = 0

	for ( i in samples:length(sa) )
	{
		subs = sa[(i+1-samples):i]
		aveSA = sum( subs, na.rm=T )/samples
		rdeltaSA = (sa[i]-aveSA)/aveSA*100
		if ( !is.na(rdeltaSA ))
		{
			if ( rdeltaSA >= mca )
			{
				freq = freq + 1
			}
		}
	}
	return(freq)
}

totalDistance <- function( x, y )
{

#total distance is in effect the same as the average velocity for the purposes of doing a single statistical test. However, for combination with other parameters where units may be important, then you may want to scale this number in your calculation. 
	
	tdistance = 0;
	for ( i in 1:(length(x)-1) )
	{
		p1 = c(x[i],y[i])
		p2 = c(x[i+1],y[i+1])
		tdistance = tdistance + distance(p1,p2);
	}
	return(tdistance);
}

timeInZone <- function (x,y,x1,y1,x2,y2)
{

#calculates the percentage of time in a
#certain boundary area.
#the most important thing to remember aboout this function is that each square is diffent, so only symmetrical areas with at least one direction of symmetry will have any meaning if you pool over all mice. Otherwise seperate the mice into groups.
	
	count = 0;
	
	minx = min(x,na.rm=T);
	maxx = max(x,na.rm=T);
	miny = min(y,na.rm=T);
	maxy = max(y,na.rm=T);
	
	lenX = as.integer(length(x));
	
	for (i in 1:length(x))
	{
		if (is.na(x[i]) || is.na(y[i]) )
		{
			#Do not count
		}
		else
		{
			if (x[i] >= minx + x1 && x[i] <= maxx - x2 &&
				y[i] >= miny + y1 && y[i] <= maxy - y2)
			{
				count=count+1;
			}
		}
	}
	
	return(count/length(x));
}

timeInZone2 <- function (x,y,x1,y1,x2,y2)
{
	count = 0;
	
	minx = min(x,na.rm=T);
	maxx = max(x,na.rm=T);
	miny = min(y,na.rm=T);
	maxy = max(y,na.rm=T);
	x1 = x >= minx +x1 & x <= maxx - x2;
	x2 = y >= miny +y1 & y <= maxy - y2;
	x3 = x1 & x2;
	count = sum(x3,na.rm=T);
	
	return(count/length(x));
}

readall <- function(k=".csv$",beginpoint=1200,dpoints=5600)
{
	ave_sa = 0;
	rearing_freq = 0;
	total_distance = 0;
	prop_in_center = 0;
	mouse_id = 0;
	edg = c_distance_from_edge;

	flist = dir(pattern=k);
	i=1;
	duration = length(flist);
	
	for (j in flist)
	{	
		#cat(i,":", j,"\n");
		k = read.table(j,sep=ftype,skip=beginpoint + meta_data_length ,na.strings="-",nrows=dpoints);

		ave_sa[i] = verticalExposure(k$V5)
		rearing_freq[i] = rearing( k$V5, r_minimum_surface_area_diff, r_samples_used)
		total_distance[i] = totalDistance( k$V3, k$V4 )
		prop_in_center[i] = timeInZone2( k$V3, k$V4, edg, edg, edg, edg )
		
		id1 = sub(mouse_number_format_remove,"",j)
		id2 = sub(".csv","",id1)
		mouse_id[i] = as.double(id2)
		
		i=i+1;
	}
	
	if (show_mouse_number)
	{
		opf_data = data.frame( mouse_id, ave_sa, rearing_freq, total_distance, prop_in_center )
	}
	else
	{
		opf_data = data.frame( ave_sa, rearing_freq, total_distance, prop_in_center )
	}
	
	write.table(opf_data, file = "opf-data.txt", sep="\t", row.names = F, quote = F)

	cat("Done.\n");
	return(0);
}

if (autorun) { readall(beginpoint = start_time, dpoints = total_time) }
